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Abstract. We study the response of the Hodgkin-Huxley neuron stimulated 
periodically by biphasic rectangular current pulses. The optimal response for 
charge-balanced input is obtained for cathodic-first pulses with an inter-phase 
gap (IPG) approximately equal 5ms. For short pulses the topology of the global 
bifurcation diagram in the period-amplitude plane is approximately invariant with 
respect to the pulse polarity and shape details. If stimuli are delivered at neuron's 
resonant frequencies the firing rate is a continuous function of pulse amplitude. 
At nonresonant frequencies the quiescent state and the firing state coexist over 
a range of amplitude values and the transition to excitability is a discontinuous 
one. There is a multimodal odd-all transition between the 2:1 and 3:1 locked-in 
states. A strong antiresonant effect is found between the states 3:1 and 4:1, 
where the modes (2 -|- 3n) : 1, n = 0,1,2,..., are entirely absent. At high 
frequencies the excitation threshold is a nonmonotonic function of the stimulus 
and the perithreshold region is bistable, with the quiescent state coexisting with 
either a regular or chaotic firing. 



PACS numbers: 87. 19.11,87. 19. In 
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1. Introduction 

Stimulation of neural fibers with a train of electric current pulses has become standard 
procedure in several areas of chnical practice. This encompasses, but is not hmited 
to, cochlear and retinal prostheses^li 12J, deep brain stimulation [J, as well as high- 
frequency conduction block of peripheral nerves [4 . Typically such devices should 
minimize charge and energy per pulse. The minimum current amplitude required for 
the neuron to spike and its dependence on stimulation frequency and pulse shape 
are among the basic properties measured in experiment. The injected charge should 
be kept at safe levels, mainly to avoid changes in the chemical content in the tissue 
near the electrode which could damage the tissue or cause the electrode corrosion [5j. 
Typically the minimum charge injected by the electrode needed by the neuron to 
spike occurs at a vanishing pulse width whereas the optimization of the energy of 
the pulse occurs at a certain nonzero width[6^i HE]. If pulses are short various, 
different waveforms have very similar efficiency 0. Differences become apparent for 
longer pulses. Most cochlear implants use charge-balanced biphasic (BP) stimulation. 
Studies of auditory nerve fibers (ANF) can be reduced [TO', HT by introducing a delay 
between standard BP pulses commonly used in cochlear implants. 

Irregular ffiing and desynchronization among a group of ANF is usually attributed 
to irregular synaptic input and physiological noise O [13]. It must be kept in 
mind however that neurons themselves are strongly nonlinear systems and are 
capable of irregular firing even in the absence of noise[Tll [T5]. O'Gorman et al. 
showed recently [T6| [17] that a dynamic instability is a plausible explanation of ffiing 
irregularities in stimulated ANF. Using the FitzHugh-Nagumo (FHN) model[18l, 
they obtained a positive Lyapunov exponent which is consistent with experiments, 
where mutual desynchronization between similarly stimulated fibers was observed. 
Electrically stimulated ANF ffie regularly at low rates and irregularly at high rates 
of stimulation. The sensitivity in this regime bears some resemblance to experimental 
results on electrically stimulated ANFfT^dD]- The FHN model is a simplification of 
the more complex Hodgkin- Huxley (HH) equations. The HH model subject to a high 
frequency monophasic pulse current also exhibits a strongly nonlinear response [14], 
including nonmonotonic dependence of the ffiing rate on pulse amplitude and a high 
coefficient of variation of interspike intervals. There is a significant slowing down of 
average response frequencies in the vicinity of the multimodal transition where the 
parity of response modes is changed. Therefore it would be useful to know if the HH 
model behaves similarly also for charge-balanced inputs. 

The present work may also help in solving problems related to the spike timing 
control in resonant neurons. The ability to control precisely the time at which neuron 
will emit a spike if a firing event occurred at some earlier time [8^ is important in some 
applications. In general the minimum time to fire is a complex issue, dependent on 
imposed constraints, initial conditions and details of the model. One of the main 
complications is the bistability of the HH model. The obtained minimum firing time 
and the minimum impulse strength may fall within a bistable regime, where a small 
change of the neuron's initial state may lead to no spike. Knowing the location of 
bistable regimes in the response diagram of a typical resonant neuron should help in 
designing spike control algorithms. Here the HH model serves as a convenient testing 
ground. 
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2. Model and results 

We will try to provide some support to earlier experimental findings using the basic 
HH model. While this model has no direct apphcation in the chnical setting we 
can improve our understanding of clinical experiments by studying the excitation 
threshold's dependence on stimulus waveform in simplified models. 

We will analyze the response of the HH neuron to periodic biphasic charge 
balanced pulses. Such stimuli were investigated experimentally [TTl and 
theoretically by several groups [TOl [TTl [IS [Mj- The requirement of charge neutrality 
is imposed to avoid tissue damage. 

We consider the model with the classic parameter set and rate constants ^5]. 

C — = —Inu — Ik - II + lapp, (1) 

where iNa, Iki II, lapp, are the sodium, potassium, leak, and external current, 
respectively. C = 1/iF/cm is the membrane capacitance. The input current is a 
periodic set of rectangular steps of height Iq and width r. The simulations are carried 
out with the time step of 0.001ms and are run for 40s, discarding the initial 4s. Fig. 
[T] shows the form of biphasic pulses used in the calculation. The emphasis of this 
work is on the charge-balanced stimulus shown in the top part of Fig. [H However the 
dependence of the threshold on IPG at different stimulation frequencies will also be 
presented for the monopolar biphasic pulses, shown in the bottom diagram of Fig. [H 
Stimuli of this form are often used in cochlear |J0, JT| |23l [24] and visual [T] [2] neural 
prostheses, as well as deep brain and muscle stimulation [5]. Electrical stimulation 
protocols have been used to learn about cortical function [26^ and vestibulo-ocular 
reflex eye movements in chinchillas [27 . High-frequency biphasic current pulses may 
also be useful in reversible peripheral nerve block that would be desirable in some 
chnical apphcations[4 . We will investigate the influence of the amplitude of the pulses 
as well as their period T^, width r and delay ri between the onsets of the two phases. 



^1 



Ti 




Figure 1. Stimulus waveform used in the calculations: charge-balanced biphasic 
form with delay (upper) and monophasic pulses with delay (lower). The delay 
between the cathodic and anodic phase is ri and r is the width of each phase. 
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Fig. [2]shows the global bifurcation diagram for a charge-balanced biphasic current 
with no delay between the cathodic and anodic part, i.e. r = ri. The firing rate 
folfii where fo is the firing frequency in response to the stimulus frequency is a 
continuous function of the current amplitude near the resonance, Ti ~ T^es — 17ms. 
Ti is the period of the stimulus. The response at antiresonant frequencies is bistable. 
Boundaries of bistable areas are found using a continuation algorithm, in which end 
values from the previous iteration are used as initial conditions for a new parameter 
set. 

The overall topology of this diagram closely resembles the result for monophasic 
monopolar pulses [28 . In the high frequency hmit the threshold is a comphcated 
function of Ti with a maximum near Ti = 5.5ms, which is approximately 1/3 of 
the resonance period Tres- Nonmonotonic behavior of hearing threshold was also 
observed in experiments in cochlear implant users [29 . Below Ti = 8ms the entire 
perithreshold region is bistable. Here the quiescent state coexists with either a locked- 
in or a chaotically firing state. Which of these solutions is obtained depends sensitively 
on initial conditions. 

The firing frequency depends nonmonotonically on Ti. Fig. |3] shows fo vs. Ti for 
/o slightly above threshold. The largest deviation from the resonant frequency occurs 
between the 2:1 and 3:1 states, near Ti = 7.5ms. This is a signature of a dynamic 
instability associated with the competition of even and odd modes (14], 
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Figure 2. Response diagram for a biphasic stimulus with no delay between the 
anodic and the cathodic part of a pulse (upper). The lower portion shows the 
high frequency limit in more detail. Borders of bistable regions are marked with 
a dotted line. The resonance near Tj = 17ms is dominated by the 3:1 state and 
higher order states. The second resonance near Ti = 34ms is occupied by modes 
higher than 1:1 (see also Ref. [28]). 

In an earlier article we have shown that the HH neuron responding to periodic 
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Figure 3. The firing frequency in the hmit of high stimulation frequency for 
a biphasic pulse at Iq = 47/iA/cm'^. The response between the mode- locked 
intervals is nonmonotonic and highly irregular. 



sequence of monophasic pulses undergoes a multimodal odd-all[TT and even-all 
transition [28^ in some regions of parameter space. Fig. [4] shows the minimum of 
the firing rate occuring between the 2:1 and 3:1 locked states for pulse amplitude 
Iq = 35/iA/cm^. In the lower part of Fig. [Uwe can clearly see that the weight of 
even modes is nonzero only for Ti above the minimum of the firing rate. For smaller 
Ti only odd modes exist. The role of the mode parity is easier to understand if we 
remember that 2Ti, where Ti ~ 8.5ms is approximately the resonant period. The 
competition of the odd and even modes in Fig. [4] accompanies the transition between 
the resonant and antiresonant regime. The vanishing of even modes below Ti ~ 8.4ms 
and a significant decrease of the firing rate are a clear signature of the antiresonant 
regime. Odd multiples of Ti = 8.4ms are associated with antiresonant stimulation 
frequencies. 
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Figure 4. The firing rate (upper) and the histogram weight of odd and even 
modes (lower) in the vicinity of the odd-all multimodal transition for a biphasic 
pulse at Iq = 35/iA/cm^. Even modes appear only above the minimum of fo/ fi- 
An interspike interval is classified as belonging to the n^^ mode if it falls between 
(n - l/2)Ti and (n + 1/2)T^. 
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Figure 5. Upper: The ratio of the output period To to the input period Ti 
between the 3:1 and 4:1 entrained states for a biphasic pulse at /q = 45.8/iA/cm^. 
Lower: Histogram of response modes at Tj = 5.55ms and Iq = 45.8/iA/cm^. The 
histogram bin size is 0.01ms. Note the absence of modes 5,8,11,14, and 17. 



Another interesting effect occurs between the 3:1 and 4:1 locked states at high 
frequencies. Bands of To vs. Ti values are shown in Fig. [5j Note the absence of all 
modes (5 + n) : 1, where n is a nonnegative integer. This is a consequence of the fact 
that the product (5 + n)r^ falls in an antiresonant regime for Ti c:^ 5.5ms. Since this is 
a qualitative effect it may provide a strong test of apphcabihty of the Hodgkin- Huxley 
model in various contexts. Quahtatively similar result should be obtained for other 
types of pulses as well. 

Fig. [H] shows the excitation threshold as a function of the delay between the 
anodic and the cathodic phases of the pulse, for four choices of the stimulus period. 
In each case the optimum ri is about 5ms. If the input period is in the antiresonance 
regime, i.e. for = 11 and 23ms, there are two thresholds because of bistability, both 
of which have minima at similar values of ri . The upper threshold where the transition 
to excitability occurs through a Hopf bifurcation, has a maximum at approximately 
Ti = 13ms. This antiresonant feature occurs at the same ri also for Ti = 34ms. Lower 
thresholds for biphasic cathodic-first stimuli were also observed in Refs. [30| [3T| |32] 
and calculated by Smit et al.[24] in an ANF model containing persistent sodium and 
slow potassium currents. The shape of the curve in Fig. [6] for Ti = 34ms is reminiscent 
of a phase response curve for the Hodgkin- Huxley model[8 . 

The minimum threshold is obtained when the stimulus is tuned to the natural 
resonance of the neuron. The maximum delay for the interaction between the cathodic 
and anodic phase to occur is approximately 18ms. The cathodic phase of the 
magnitude shown in the bottom diagram of Fig. |6] induces subthreshold oscillations 
with a resonant period. The spike-triggered average of stimulus I{t) for isolated spikes 
in the HH model driven by Gaussian random noise current with a short correlation 
time also indicates the preference of the neuron for approximately 5ms separation 
between the negative and positive parts of the current^. This is a reflection of the 
HH neuron's internal dynamics in which the sodium channel remain open for about 
5ms from the stimulus onset. The bifurcation diagram with optimal IPG, ri ~ 5ms, 
closely resembles Fig. [H The only significant difference is a shift to lower values of Iq. 

The threshold for biphasic monopolar pulses is shown in Fig. [8l For Ti < Tres 
the most regular input signal with equally spaced pulses is also the least hkely to ehcit 
spiking. The meaning of parameters is as in Fig. [H When the regularity of the signal 
is distorted the threshold decreases. 
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Figure 6. Excitation edge for different interpulse separation as a function of 
time difference between the cathodic and anodic phase. Prom top to bottom: 
Tj = 11, 17, 23, and 34ms, for stimulation by rectangular current pulses of width 
T = 0.6ms and height Iq. For different pulse frequencies the minimum occurs at 
ri 5ms separation between the cathodic and the anodic phase. The threshold 
is bistable for = 11ms and 23ms. 
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Figure 7. Response diagram for stimulus with an optimal inter-phase gap. 
Bistable regions are marked b. The resonance area near Tj = 17ms is dominated 
by the 3:1 state and higher order states. The second resonance near Ti = 34ms is 
occupied by modes higher than 1:1. 



Response of the Hodgkin- Huxley neuron to a periodic sequence of biphasic pulses 8 




5 10 
15 I 1 1 .— 




5 ^ ' ' ' ^ 

5 10 15 20 
15 I . . . ^ ^ r- 



10 - 




5 ' ^ ' ' ' ' ' ' 

5 10 15 20 25 30 

(ms) 

Figure 8. Excitation edge for a double pulse as a function separation between 
the leading and the trailing phase. The diagrams from top to bottom are drawn 
for Ti = 11, 17, 23, and 34ms, respectively. 

Comparing Fig. [7] to Fig. [2 we can conclude that the topology of the response 
diagram does not depend on the size of IPG for ri < 5ms. The same is hkely to be true 
for all types of periodic short pulses irrespective of details of their time dependence. 
This is an illustration of the well known fact that the HH neuron has a charge threshold 
property. When a relatively strong current is delivered in a short time, the voltage 
change is mainly determined by the capacitive current [34l [35l [36] . 

In neural electrical stimulation the main safety factors are the amount of charge 
transferred in a single pulse, Qth, and energy per spike, Eth, needed to evoke spiking, 

Qth'- I Ith{t)dt (2) 
Eth- r I?Hit)dt, (3) 

^0 

where Ith is the threshold current amplitude. Fig. [9] shows the threshold charge at 
Ti = 17ms as a function of the phase width r for ri = 5ms. The charge-duration 
curve in Fig. [9]imphes the use of short pulses in stimulation protocols if the injected 
charge is to be minimized. In practice, pulses lasting tens of microseconds approach 
the minimum charge condition sufficiently well and are often a reasonable solution in 
the design of neural prostheses. During this relatively short time one may be able 
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Figure 9. Threshold charge as a function of the r for Tj = 17ms and ri = 5ms. 
The minimum occurs at r = which means that the most effective pulses are 
those with vanishing width. 

to avoid Faradaic reactions that would occur at higher levels of total charge with 
longer pulses. Similar conclusion was reached in a study of Sahin and Tie|7J who 
investigated effects of non-rectangular waveforms on threshold both experimentally 
and theoretically within a simple local model of a mammalian nerve [33 [38]. Their 
model is similar to the HH model expcept for the absence of the potassium channel. 
However if the main hmiting factor is the energy delivered per pulse, then finite r 
is preferable, see Fig. [TOl Qualitatively similar result was obtained for different 
waveforms in Ref. 




T (ms) 

Figure 10. Threshold energy at the resonance frequency, Ti = 17ms, and 
ri = 5ms, as a function of the time separation between the leading edge of the 
cathodic phase and the leading edge of the anodic phase. The minimum occurs 
near r = 4ms. 



3. Conclusions 

Studies of threshold behavior of resonant neurons may help designing future cochlear 
implants. Nonlinear instability present in the HH model may provide important clues 
to understanding ANF response to high-frequency stimulation. O' Gorman et al. used 
the FHN model to explain irregular firing of auditory nerve fibers and the lack of 
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synchronization of response to external drive [TB]. This instabihty may be related to 
the recently discovered parity transitions of response modes [HI [28]. It is important to 
realize that a resonant cell such as the HH neuron fires regularly in resonant regimes 
and irregularly in antiresonant regimes. In our work on the HH model stimulated by 
monopolar monophasic pulses we found the odd-all multimodal transition ||T4| [15], in 
which the vanishing of even response modes is accompanied by a local minimum of the 
firing rate. This effect occurs between the 2:1 and 3:1 states and is most pronounced 
near excitation threshold. Here we found the same phenomenon in the HH neuron 
driven by charge-balanced biphasic current pulses. We also discovered a new effect: 
the vanishing of the modes 2, 5, 8, . . . in the regime of irregular firing between 3:1 and 
4:1 states. These antiresonant effects are quite general and do not depend on details 
of pulse shape or polarity of the current pulse. 

The firing threshold of the HH model is a nonmonotonic function of frequency 
with the minimum at the main resonance frequency. A local minimum exists also 
at very high stimulation rate. At the resonant frequencies the solution is always 
unique and the firing rate is a continuous function of the pulse amplitude. At non- 
resonant frequencies the HH neuron is bistable, i.e. the firing state coexists with the 
excitable state, and the transition to excitabihty is discontinuous. Charge-balancing 
has no effect on the topology of the global bifurcation diagram in the period-amphtude 
plane. Fig. [2]closely resembles results obtained for monophasic stimulation [28]. These 
properties are quite general, provided the width of the current pulses is sufficiently 
small, which implies ri < 5ms for the classic HH parameter set, and should hold for 
other, non-rectangular pulse shapes. 

Results in Figs. [GlandfTOlindicate that IPG giving the lowest excitability threshold 
of the HH membrane is approximately 4 — 5ms. The minimum threshold approaches 
that of the excitatory monophasic current pulse and is obtained for ri ~ 5ms. 
The decrease of threshold associated with the increase of IPG is consistent with 
experimental measurements [TO] . Similar result was obtained for biphasic pulses with 
long IPG in a phenomenological model of the electrically stimulated human ANF[23j. 
Our results are also consistent with the work of Carlyon et al.[ll^ who studied human 
behavioral thresholds for trains of biphasic pulses apphed to a single channel of 
cochlear implants as a function of IPG. The threshold decreased for IPGs up to several 
milliseconds when the phases of the pulse were of opposite polarity. The IPG of 
about 4ms also approximately minimizes the energy delivered per pulse. The optimal 
stimulation approach may include a Gaussian noise superposed on the input signal[39]. 
Small amounts of noise were shown to improve detection of faint signals, improve 
temporal resolution and regularize strongly irregular firing in regimes of dynamic 
instabilities. 
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